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In this paper we present two new numerical methods for studying thermodynamic quantities 
of integrable models. As an example of the effectiveness of these two approaches, results from 
numerical solutions of all sets of Bethe ansatz equations, for small Heisenberg chains, and Monte 
Carlo simulations in quasi-momentum space, for a relatively larger chains, are presented. Our results 
agree with those obtained by thermodynamics Bethe ansatz (TBA) and Quantum Transfer Matrix 
(QTM). 

PACS numbers: 75.10.Jm, 75.40.-s 



I. INTRODUCTION 

The study of exactly solvable models is a very impor- 
tant field in condensed matter physics, which began with 
Bethe's solution of the isotropic Heisenberg chain Q. In 
general, the Bethe ansatz (BA) solution of a model has 
several drawbacks: it has a complex mathematical struc- 
ture; the excitations are not immediately available; and 
most important, it does not give explicit results even for 
the thermodynamic quantities of the system. It was only 
when Yang and Yang Q presented a strategy to study 
the thermodynamics of BA solvable systems that the 
temperature dependence of quantities such as the spe- 
cific heat and the magnetic susceptibility become avail- 
able. The method is now designated as thermodynamic 
Bethe ansatz and has known many developments in the 
last thirty years ||. Additionally, correlation functions, 
such as the conductivity, can not be obtained from the 
BA equations alone, and combination of BA results with 
other methods is required for their calculation JI]. 

The BA ansatz method has been applied to Bose, 
Fermi jj| ^, || , and spin systems fl fjj . It is a general fea- 
ture of the BA solution, first proved by Yang and YangQ 
for the Bose case, that a given eigenstate of the model is 
characterized by a unique set of quantum numbers {Ij}- 
Further, it also can be shown that all configurations of 
these quantum numbers Ij exhaust the Hilbert space of a 
given model. Since the energy eigenvalues are functions 
of the above quantum numbers, instead of using TBA 
and quantum Monte Carlo approaches, we can study BA 
solvable models in quantum number space by classical 
Monte Carlo. Furthermore, for a small system (but for 
larger systems than those available to exact diagonaliza- 



tion methods), it is possible to solve the BA equations 
for all eigenvalues. Therefore, the expectation value of an 
Hermitian operator in thermal equilibrium can be com- 
puted. 



In this paper, we shall introduce two numerical ap- 
proaches for computing thermodynamic quantities of 
Bethe ansatz solvable models. The methods are ilus- 
trated with the ID isotropic Heisenberg model, since 
this model is well studied in the literature. Further- 
more, the study of the Heisenberg model is it self rele- 
vant, since this system predicts many properties of quasi- 
one-dimensional materials [|l0| O, [l2|. This model has 
been investigated by many kinds of methods. For exam- 
ple, the low temperature behaviors are quite well under- 
stood by a combination of Bethe ansatz Jl3[ and con- 
formal field theory ju| |l^|. A strong logarithm singu- 
larity in the susceptibility at low temperature was first 
found by the Bethe ansatz calculation of quantum trans- 
fer matrix ]l6| and then verified experimentally Jl0[ pd| . 
The thermodynamics of the model has been studied by 
TBA 




L9 20L El 221 as well as by QTM 
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The paper is organized as follows. In Sec. [n[ we first 
briefly review the BA solution of the isotropic Heisen- 
berg model. In Sec. |n| and Sec. IV, we introduce the 
basic idea of numerical Bethe ansatz (MBA) and Monte 
Carlo Bethe ansatz (MCBA). In Sec.[y| we check the ef- 
fectiveness of these two methods computing the specific 
heat and the magnetic susceptibility in the absence of 
an external magnetic field and compare our results with 
those obtained from the TBA. We than use our methods 
to study the two quantities above in the presence of an 
external magnetic field. A brief summary is given in Sec. 
VI. 
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II. ISOTROPIC HEISENBERG MODEL 

Now let us first review the Bethe ansatz solution of the 
ID Heisenberg chain, which can be found in the book of 
Takahashi |3| . The Hamiltonian of the isotropic Heisen- 
berg model is 



JY 



(1) 



1=1 



where N is the number of sites, S?,Sf,S? are spin 



1/2 operators at site I and J 



J i ' u i 
-1,1 represent anti- 



ferromagnetic and ferromagnetic cases, respectively. The 
solution with periodic boundary condition Sat+i = Si 
using the string hypothesis takes the form 



N9(x™/ri) = 2nl^ + ]T 9 nm (a 



Here 6{x) — 2 tan and 



•a#). (2) 



e nro (x) = 9 ( , - , )+20[ , )+••• + 20 I — — ? ) + ( — f— ) For : n ^ m 



|n — m|/ \Jn — m| + 2/ \n + m — 2J \7i + m 

4/ l2n-2/ V2n 



20 (T) + 26» f ^) +... + 20 ( — ^— J +0(^-1 For:n = m (3) 



and x™ is the real part of the n-string which is designated 
by 



x^=x 7 + i(n+l-2j), j = l, 



, ,71 



I" is the quantum number of 7th n-string (note that n 
and 7 are indices). We denote the number of the n-string 
by a n , thus 7 = 1, . . . , a n and the string configuration 
{a} satisfy 

ai + 2a 2 H 1- (M - 1)«m-i + Ma M = M, (4) 

where M is the number of down spins. The quantum 
number of n-string 7™ is an integer (half-odd integer) if 
N — a n is odd (even) and satisfy 



|J™| < (JV — 1 



E 

m— 1 



(5) 



where t nm = 2 min(n, m) — S nm . For a given set of {12}, 
Eq. (^|) can be solved numerically and the energy is given 
by 



E{I2} = -NJ/A + Y, 



2Jn 



n,j 



(x™) 2 +n 2 



(6) 



which represents the energy of the lowest weight state in 
SU(2) irreducible space designated by S = N/2—M, S z = 
S, S — 1, . . . , — S. In the presence of an external field h a 
Zeeman term is added to Eq. @. Hence the total energy 
of a given quantum number configuration is given by 

E = E{I^}-hM, (7) 

where M. = 2S Z is the magnetization of the state 



III. NUMERICAL BETHE ANSATZ 

In statistical mechanics, the expectation value of an 
Hermitian operator Q in thermal equilibrium is given by 



(8) 



where Z is known as partition function, defined as 



[3 is inverse temperature, and represents sum over 
all possible eigenstates of the Hamiltonian. It turns out 
that the variation of Z with respect to temperature or 
any other external parameters affecting the system can 
virtually tell us everything we might want to know about 
the macroscopic behavior of the system. For example, the 
internal energy is given by 



(10) 



From Eq. (|9|), it is easy to see that the internal energy can 
also be written in terms of a derivative of the partition 
function: 



U 



1 dZ _ d\nZ 
Z~df3 ~ 0/T' 



(11) 



The specific heat is given by the derivative of the internal 
energy: 



_ dU , a2 dU , a2 d 2 \nZ 



(12) 
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where ks is the Boltzmann constant which is set to unit 
hereafter. 

Our aim is to combine the idea of statistical mechanics 
mentioned above with the numerical solution of the BA 
equations. The main idea of the numerical Bethe ansatz 
method we introduce here is, first, to compute all eigen- 
values of a BA solvable model from its corresponding BA 
equations. Then to compute the expectation value of the 
Hermitian operators, representing the physical observ- 
ables we are interested in, by averaging those operators 
over all states of the system, weighting each state with 
its own Boltzmann weight. 

It has been shown || that the Hilbert space of the 
isotropic Heisenberg model is complete under the string 
classification. Here we want to show how to travel 
through all C^ 2 states in quantum number space and 
illustrate it by considering a system of 6 sites. 

For the case of M down spins, the first task is how to 
obtain all string configurations fulfilling the restriction 
(|J). We adopt a time-like number u aM ■ Qjj-i : ••• : 
a 2 '■ o-i" , where the magnitude a n measures from to 
[M /n] (here [a;] returns the truncated integer value of x) , 
just like hour and minute in "HH:MM" measure from 
to 23 and to 59 respectively. If we increase the number 
"ajvf : • ■ • : a{" by adding 1 to the first digit a\ step by 
step, we can travel through all possible values. Among 
all these numerical values, only those whose digits satisfy 
the condition (Q) are what we need. Then all string con- 
figurations can be found by this procedure. Of course, 
these operations are realized in a computer. In order to 
make the method clear, let us consider a problem of 6 
sites. 

M=0: It is easy to have the state with all spins up, i. 
e. M = 0, which has energy E = —JN/A. 

M=l: In this case, we only have one string configu- 
ration ai = 1 and one quantum number —2 < Ji < 2, 
thus there are 5 states. Each of them is represented by 
one quantum number in the interval [—2, 2]. We can get 
all possible quantum number configurations from the fol- 
lowing figure, 

o o • o o 

where the dot is the occupied quantum number, and the 
open circles represent other possible quantum numbers. 
Then the BA equation is just 

6 tan -1 x\ — nli, (13) 

which has a simple solution x\ = tan(7r/i/6). 

M=2: Here the string configuration is characterized 
by {a\,a 2 }. We construct a number "a 2 : a\ \ in which 
the maximum value of ot\ is 2 ([2/1] = 2), and a 2 1 
([2/2] = 1). Increasing ol\ step by step we generate all 
possible configurations of the "«2 : «i" number, ranging 
from 0:0 to 1:2. Among all the generated configurations, 
we are only interested in those satisfying the condition 
a,\ + 2a 2 = 2. The first case is a.\ = 2, a 2 — 0, in 
which the quantum number satisfy —1.5 < l\,l\ < 1.5, 



TABLE I: All quantum number configurations for M — 2. 



a\ — 2, 012 = 


n 
n 


-1.5 
-0.5 


-1.5 
0.5 


-1.5 
1.5 


-0.5 
0.5 


-0.5 
1.5 


0.5 
1.5 


CKl = 0, OLI = 1 


n 


-1 





1 









the second one is ot\ = 0, ot 2 = 1, in which the number 
satisfy — 1 < if < 1. They can be characterized by 

o • • o 

and 

o J o 

respectively, where J denotes the occupation for a quan- 
tum number of 2-string. In Table |, we list all quantum 
number configurations for M — 2. The BA equations for 
these two cases are 

6 tan" 1 x\ = nil + tan_1 ^ ~ 2 — i 

6 tan" 1 x\ = nil + tan_1 ^ ~ 2 ^ > (I 4 ) 

and 

6tan" 1 (^/2) = nlf, (15) 

respectively. 

M=3: In this case the string configuration is char- 
acterized by ct2, a^}. In the same way as we did 
above, we construct a number "«3 : a 2 ■ c*i \ the max- 
imum value for each digit from left to right is 1, 1, 3 
respectively. Then we have 3 string configurations with 
the condition a\ + 2a 2 + — 3, which correspond to 
the following sequences 

a : • • • 



6: o • o 



where a, 6, c have 1, 3, 1 states respectively, § denotes 
the site for 3-string. And in Table [n], we list all quantum 
number configurations for M = 3, whose BA equations 
are 



6 tan - 




= nl\ + 
+ tan 


tan 1 (x\ 
-\x\-x\ 


— X 
), 


6tan~ 


x 2 


= + 
+ tan 


tan -1 (^2 

(x 2 — x. 


— X 

), 


6 tan" 


1 T 1 

^3 


= + 
+ tan 


tan -1 {x\ 
1 (x\ — X 


— X 

)• 



4 



TABLE II: All quantum number configurations for M = 3 



O'l 


— 3, Ct2 


= 0, a 3 


= 


n 


-1.0 














I\ 

















II 


1.0 






«r 


= 0, a 2 


= 0, a 3 


= 1 


I\ 


-1.0 





1.0 










n 













= 1, «2 


= 1, a 3 


= 


it 










the Monte Carlo calculation, which follows three steps. 
Let us assume the present state is /i with a correspond- 
ing Mfj, - the number of down spins in state /i. From the 
state n any other state v with M v , in number of C^y 2 — 1, 
can be obtained. 

step one: first we choose M u> knowing that the num- 
ber of states with M v spins down is CjJ£ — Cm„-i> tnus 
the probability of selecting M v is (C^ — C^_ x )/C^, 2 . 

step two: selected M v , all possible string configura- 
tions given M v are determined from of Eq. (|]) which 
satisfy § 



6 tan 1 x\ 



6 tan 1 x\ — 



nil + tan _1 (a;} - x\) 
+ tan- 1 ((x 1 -x\)/Z), 
nil + tan -1 ^ - x\) 
+ tan- 1 ((^?-^i)/3), 



and 



respectively. 



6 tan 1 . 



nil 



(17) 



(18) 



As a result we have totally Cf = 20 distinct configu- 
rations of quantum number whose Hilbert space is com- 
plete. 

Then, we compute the eigenvalue for a given quantum 
number configuration {/"} by solving the BA equations 
numericaly. For the Heisenberg chain, the BA equations 
can be solved by iteration, for other models, such as the 
Hubbard model, the BA equations can be solved by a 
gradient method. 



IV. MONTE CARLO BETHE ANSATZ 

For a system of N sites, there are C^, 2 quantum num- 
ber configurations. This number increases exponentially 
with the size of the system, so it is impossible to calcu- 
late all eigenvalues for a large system, such as N > 40, 
under present computer capacity. This restriction can be 
overcome by a Monte Carlo method. There are many 
Monte Carlo methods available, and we introduce below 
a new method that we call Monte Carlo Bethe ansatz. 
This method is a classical Monte Carlo strategy applied 
to a quantum problem. The basic idea behind the MCBA 
method is to simulate the random thermal fluctuation of 
the system from state to state in quantum number space 
of the BA solution. This method is not limited by the 
sign problem, that may show up in the usual quantum 
Monte Carlo methods. 

Since the energy eigenvalues are a function of both M 
and of the quantum numbers /" we can follow a classi- 
cal Monte Carlo strategy, by sampling the configuration 
space of M and {12}. We now explain how to implement 



E 

aiH hMa fl 



D({a n }) = (C£ - C 



M-l) > 



(19) 



with D({a„}) is the number of states, characterized by 
the set of quantum numbers associated with the 

string configuration {«„}, and reads 



D({a n }) = l[C aj 
»=i 



(20) 



So, in step two, we select a string configuration with the 
probability D({a n })/(C» - C&_ x ). 

step three: having determined the string configura- 
tion, we then select at random a quantum number con- 
figuration, which is the state v we want, for the given 
string configuration. From the partition function the 
probability density for a state fj, is 



p M = (N — 2M„ + l)e 



(21) 



where the degenerancy of state [i was taken into account. 
The detailed balance condition tells us the transition 
probability should satisfy 



Pu = (N - 2M U + 1) g(Ev _ Ea) 
p„ (N - 2M^ + 1) 



(22) 



Hence it is possible to use the Metropolis algorithm for 
the acceptance ratio to accept or reject the state fi ac- 
cording to 



A(p -»!/) = 



-2A/,„ + l) 0(E„-E,A 



-2M M + 1)' 



< 1 



pI " \ (23) 
otherwise. 



The MCBA algorithm is complete and the three basic 
steps are repeated a number of times. After an initial 
equilibration time, the expectation values can be then es- 
timated as an arithmetic mean over the repeated Markov 
chain 



(Q) 



(24) 



{/-} 
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FIG. 1: The specific heat of 24 sites anti- ferromagnetic and 
ferromagnetic XXX model (points) and the same quantities 
obtained by TBA (lines). 




■rial 



FIG. 2: The susceptibility of 24 sites anti-ferromagnetic (left) 
and ferromagnetic (right) XXX model (points) and the same 
quantities obtained by TBA (lines). 



z i i i i i i i ; 




T/IJI T/IJI 

FIG. 3: The specific heat and susceptibility for a 60 sites 
chain, computed with the MCBA method, is compared with 
the TBA results (lines) both for the anti-ferromagnetic (cir- 
cles) and ferromagnetic (squares) cases. 



TABLE III: Specific heat and susceptibility of ferromag- 
netic XXX model obtained by thermodynamic Bethe ansatz 
(TBA), numerical Bethe ansatz (NBA) solution of 24 sites, 
and Bethe ansatz based Monte Carlo (MCBA) approach for 
60 sites system. 





T/J 


TBA 


NBA 


MCBA 


Cv 


0.5 


0.129178 


0.129197 


0.1261 ± 0.0019 




0.6 


0.121671 


0.121688 


0.1220 ± 0.0012 




0.7 


0.112363 


0.112379 


0.1132 ± 0.0009 




0.8 


0.102496 


0.102511 


0.1019 ± 0.0003 




0.9 


0.092861 


0.092875 


0.0931 ± 0.0002 




1.0 


0.083875 


0.083887 


0.0840 ± 0.00017 




1.5 


0.051111 


0.051117 


0.05112 ± 0.00004 




2.0 


0.033256 


0.033260 


0.03326 ± 0.00002 




2.5 


0.023088 


0.023090 


0.02308 ± 0.00001 




3.0 


0.016876 


0.016878 


0.01687 


X 


0.5 


3.7378 


3.742446 


3.669 ± 0.027 




0.6 


2.90686 


2.909917 


2.8914 ± 0.0139 




0.7 


2.35856 


2.360709 


2.3458 ± 0.0039 




0.8 


1.97323 


1.974817 


1.9613 ± 0.0044 




0.9 


1.68953 


1.690744 


1.6902 ± 0.0028 




1.0 


1.473032 


1.473986 


1.47509 ± 0.00180 




1.5 


0.882613 


0.882996 


0.88316 ± 0.00016 




2.0 


0.622855 


0.623056 


0.62299 ± 0.00010 




2.5 


0.479082 


0.479205 


0.47908 ± 0.00007 




3.0 


0.388433 


0.388516 


0.38848 ± 0.00004 



V. SPECIFIC HEAT AND SUSCEPTIBILITY 

In order to check the validity of our approaches, we 
apply these two methods to the study of the specific 
heat and of the magnetic susceptibility of the anti- 
ferromagnetic and ferromagnetic Heisenberg models. 

For the present model, however, because of the degen- 
eracy in each set of quantum number configuration, Eq. 
(||) should be revised according to the property of the 
operator. For example, the internal energy and magneti- 
zation are 

{E) = ±j2( N - 2M » + 1 ^ e ~ PE ^ 

< M > = ^E E 2Mle-f> E ». (25) 

where Z = J2^( N ~ 2M v + l)e -/3£ ^. From thermody- 
namics it is easy to have the expression for specific heat 
and magnetic susceptibility per site 

C=^{{E*)~{E)% 

X = ^((M 2 )-(M) 2 ) (26) 

We apply NBA to 24-site system and MCBA to 60-site 
system, respectively. The later has C|q different quan- 




FIG. 4: The specific heat of the anti-ferromagnetic Heisenberg model for different values of the external field: (a) 
h — 0, 0.2, 0.3, . . . , 1.0; (b) h — 1.1, 1.2, . . . , 2.0; (c) specific heat as a function of h for different temperature T = 0.1, 0.2, . . . , 4.0. 
In panel (a) of fit to the law C —oc T 1//2 , for h=1.0, is given at low temperatures. 




FIG. 5: The magnetic susceptibility of the anti-ferromagnetic Heisenberg model for different values of the external field: 
(a) h — 0, 0.4, 0.5, . . . , 0.9; (b) h = 1.0, 1.2, . . . , 2.0; (c) magnetic susceptibility as a function of h for different temperature 
T = 0.1, 0.2,..., 4.0. 



turn number configurations, hence it is impossible to cal- 
culate all the eigenvalues of the system. 

In Figs. [I] and ||, we show the specific heat and the 
magnetic susceptibility, for a 24-site system, obtained 
from NBA and compare our results with those obtained 
from TBA. It is clear that the two results match. In Fig 
|§|, we show the specific heat and the magnetic suscepti- 
bility, for a 60-site system obtained from MCBA together 
with the results from TBA. They both agree to each other 
except at low temperature. In Table [v|, we compare, for 
the ferromagnetic case, the two methods we introduced 
here with TBA, giving the explicit numerical values. It 
is clear that our methods work very well for the present 
model. Hence our conclusion is that for small systems, 
such as N < 38, and for the Heisenberg chain, it is possi- 
ble to compute all eigenvalues and to obtain all possible 



thermodynamic quantities of interest by using Eq. (g). 
For temperatures larger than the finite size energy gap 
our results agree with TBA results exactly. For larger 
systems, however, results can still be obtained by using 
the MCBA method. 

Now we study the thermodynamics of the model in 
the presence of a magnetic field by NBA, which has also 
been studied by Klumper |2(|. In Fig. |] and Fig. [jj, 
the results for the specific heat and the magnetic sus- 
ceptibility of the anti-ferromagnetic case are shown for 
various magnetic fields. It is clear from these two figures 
that there are two different behaviors at low tempera- 
ture, separated by the saturation field h c — 1.0 at the 
ground state. In order to understand better this behav- 
ior of the antiferromagnetic case, let us use the mapping 
between the Heisenberg model and the spinless fermion 




FIG. 6: The specific heat(a) and the magnetic susceptibility(b) ferromagnetic Heisenberg model as a function of temperature for 
different values of external field h — 0.0, 0.1, 0.2, . . . , 2.0, obtained by NBA. (c). The susceptibility of ferromagnetic Heisenberg 
model for h = 0.0, 0.1, 0.2, obtain by NBA and MCBA respectively. 



model. This mapping is achieved by the Jordan- Wigner 
transformation [ p9| , and Hamiltonian (|l]) can be written 
as 

j N 

i=i 

N i i 

+^E^-2 )( " i+1 ~2 ) ' (27) 
i=i 

where the spinless fermion operators / ; , fi obey the usual 
anti-commutation relation, m is the usual local number 
operator. When h < h c , the system is not fully polarized, 
that is Yli=i n i > 0j hence we always have two Fermi 
points ±fcp at the ground state. The dispersion rela- 
tion of low-lying excitations is dominated by the linear- 
k dependence, hence we still have the Fermi-liquid like 
specific heat: C oc T at low temperatures. If h > h c , 
however, and from the point view of spinless fermions, 
we have n i = 0, and the dispersion relation be- 

comes fc 2 , because of the cosfc dispersion-relation for the 
fermions in the lattice. Hence, the specific heat mani- 
fest a T 1 / 2 behavior at sufficiently low temperature for 
h = h c , which can be seen in Fig. |] (open circles). More- 
over, the magnetic susceptibility presents a strong peak 
for h = h c , when T — > [see Fig. § panel (c)]. This 
strong magnetic response is associated with a change in 
the nature of the elementary excitations when the line 
h = h c is crossed at zero temperature. Indeed, at T = 
and h c = 1.0, the system manifest infinite susceptibil- 
ity, as can be seen from Fig. |J panel (b). We attribute 
it due to the degeneracy between the state of [N — 1, 1] 
and [N] , and a small magnetic field can fully polarize the 
system. The phase with h > h c share anti-ferromagnetic- 
like behavior [Fig. ||, panel (b)], while for h < h c , the 
susceptibility shows a logarithm singularity p6| . 

For the ferromagnetic case the specific heat and the 
magnetic susceptibility are ploted in Fig. pi for different 



values of the magnetic field. As is known, if ft. = the 
ground state of the ferromagnetic case is highly degen- 
erate with S = N/2, S z = -S, -S + !,■■■ ,S and a very 
small h can fully polarized the system. So it is easy to 
understand why zero temperature susceptibility is infi- 
nite. After it is magnetized (in the presence of small h), 
however, the susceptibility should be zero. This behav- 
ior is seen in Fig. ^], panel (b). We also show, in Fig. 
^, panel (c), the susceptibility obtained by MCBA. Both 
the results of the two methods agree with each other per- 
fectly. 



VI. SUMMARY 

In summary, we presented two numerical approaches 
to discuss the thermodynamics of Bethe ansatz solvable 
models. The first one is the numerical Bethe ansatz 
which works very well for a small system. We think it is 
possible to obtain all eigenvalues of a system up to size 
L = 38, for the Heisenberg model. For a relatively larger 
system, we also find that the Monte Carlo simulation 
in quasi-momcntum space works well in the moderate 
and high temperature regions. At low temperatures the 
present selection method is not excellent, and a better 
one is required. The discovery of such a method is a 
challenging and interesting research problem. 

There are many physical quantities of interest at finite 
temperature which are still not well understood, such as 
spin stiffness of XXZ model, important to understand the 
transport properties, because of the complex form of the 
thermodynamic equations. Our methods provide a new 
route to compute all these quantities directly from the 
Bethe ansatz equations. 
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Cheung Kong projects of China Education Ministry and 
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